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We present an ah initio analysis for the ground-state properties of a correlated organic com- 
pound /t:-(BEDT-TTF)2Cu(NCS)2. First, we derive an effective two-dimensional low-energy 
model from first principles, having short-ranged transfers and short-ranged Coulomb and ex- 
change interactions. Then, we perform many-variable variational Monte Carlo calculations for 
this model and draw a ground-state phase diagram as functions of scaling parameters for the 
onsite and off-site interactions. The phase diagram consists of three phases; a paramagnetic 
metallic phase, an antiferromagnetic (Mott) insulating phase, and a charge-ordered phase. 
In the phase diagram, the parameters for the real compound are close to the first-order Mott 
transition, being consistent with experiments. We show that the off-site Coulomb and exchange 
interactions affect the phase boundary; (i) they appreciably stabilize the metallic state against 
the Mott insulating phase and (ii) enhance charge fluctuations in a wide parameter region in 
the metallic phase. We observe arc-like structure in Fermi surface around the region where the 
charge fluctuations are enhanced. Possible relevance of the charge fluctuations to the experi- 
mentally observed dielectric anomaly in the ^c-BEDT-TTF family compounds is also pointed 
out. 

KEYWORDS: organic conductors, first principles, Hubbard-type low-energy model, variational Monte Carlo 
method, Mott transition 



1. Introduction 

In organic conductors, we find diverse properties and 
a plenty of phases including normal metals, supercon- 
ductors and various types of insulators with antiferro- 
magnetic, charge or (spin) Peierls orders as well as spin- 
liquid- type nonmagnetic Mott insulators.^ Although the 
unit cells of these conductors contain many atoms con- 
stituting the molecules with complicated crystal struc- 
tures, band structures near the Fermi level are in most 
cases simple with a small number of bands isolated from 
other bands located away from the Fermi level. These 
isolated bands originate from molecular orbitals [lowest 
unoccupied molecular orbitals (LUMO) and highest oc- 
cupied molecular orbitals (HOMO)]. Because the number 
of bands near the Fermi level is small, screening of the 
electron-electron Coulomb interaction is poor. In addi- 
tion and more importantly, large lattice constants make 
the overlap of the neighboring molecular orbitals small, 
leading to a large ratio of the screened electron interac- 
tion to the kinetic energy. This is the reason why electron 
correlations are in general strong in the organic conduc- 
tors. Because of the complex unit cell and prominent 
strong correlation effects, the ah initio calculation and 
clarification of mechanisms of material properties in the 
organic conductors remain as big challenges. 

Among all, a family of compounds, (ET)2A' with 
a number of choices of anions X alternatingly 
stacked with BEDT-TTF molecules [where BEDT-TTF 
is bis(ethylenedithio)-tetrathiafulvalene, abbreviated as 
ET] , offers a variety of prototypical behaviors of strongly 
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correlated electron systems with two-dimensional (2D) 
anisotropics.^ In particular, the /^-type compounds char- 
acterized by dimerization of ET molecules have served 
to discoveries of unconventional quantum phases and un- 
explored concepts at the forefront of condensed matter 
physics. 

Unconventional superconductivity is found in some of 
these compounds that show metallic behaviors. Namely, 
the compounds with the anions X=Cu[N(CN)2]Br 
(ref. 2) and X=Cu(NCS)2 (refs. 3, 4) abbreviated as 
n-^T and /t:-NCS hereafter, respectively, show supercon- 
ducting transitions at Tc ~10-13K. Under pressure, the 
compound with X =Cu2(CN)3 referred to as /^-CN also 
shows superconductivity below 2.8 K.^ However, the 
driving mechanism of the superconductivity is not fully 
understood yet. 

An unconventional nonmagnetic Mott-insulating 
phase found near the Mott transition for a>:-CN is 
another example of such a discovery. In contrast to a 
naive expectation for a magnetic order at sufficiently 
low temperatures, no apparent magnetic orders are 
observed even at a prominently low temperature T=0.03 
K that is four orders of magnitude lower than the 
antiferromagnetic spin-exchange interaction J~250 K.^ 

The emergence of the quantum spin liquid near the 
Mott transition on two types of 2D Hubbard models 
with geometrical frustration effects has already been pre- 
dicted in earlier numerical studies with the help of es- 
sentially exact algorithm of the path integral renormal- 
ization group (PIRG).^^'^^ The ground state does not 
seem to break symmetries so far proposed in literatures 
because of the geometrical frustration, while the full un- 
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derstanding of the spin liquid phase remains a chahenge. 
Most of later numerical^^ and theoretical^^ studies have 
also been performed for a simplified single-band 2D Hub- 
bard model based on an empirical estimate of parameters 
by following extended Hiickel calculations.^^' We cer- 
tainly need a more realistic and ab initio description of 
hz-ET compounds to estabhsh the existence of such truly 
as-yet-unestabhshed states and to get insights into pos- 
sible fundamentally new ideas. 

Another seminal finding achieved in this family is the 
unconventional character of the Mott transition found 
for X=Cu[N(CN)2]Cl (abbreviated as K.-CI) under pres- 
sure. The novel universality class observed by the resis- 
tivity at this Mott transition is in good agreement with 
the marginal quantum criticality at the meeting point 
of the symmetry breaking and topological change. ^^^^ 
Its significance to physics and theoretical concept of the 
quantum criticality calls further experimental test and 
critical examinations based on the comparison with the 
realistic and first-principles grounds. 

In spite of these innovative ideas and findings, the ab 
initio studies are so far few^^ and most of the studies 
were performed using empirical models inferred from the 
Hiickel studies. In addition, the strong electron correla- 
tion in the organic conductors hardly justifies its naive 
applications of the standard density functional theory 
(DFT). 

To overcome the limitation and deficiency of the con- 
ventional DFT approach, a hybrid first-principles frame- 
work, which we call Multi-scale Ab initio scheme for 
Correlated Electrons (MACE) has been developed and 
applied to a number of strongly correlated materials. 
Our general framework consists of (1) ab initio calcu- 
lations of the global electronic band structure either 
by the density functional theory or by other many- 
body theory such as the GW method and (2) a sub- 
sequent downfolding procedure by elimination of degrees 
of freedom far away from the Fermi level, which gener- 
ates low-energy effective models. ^^^^ It is followed by 
(3) the procedure to solve the low-energy models by 
more reliable low-energy solvers such as PIRG, dynami- 
cal mean- field theory (DMFT), 2^ and the many-variable 
variational Monte Carlo (mVMC) methods. The accu- 
racy of this three-stage scheme has critically been tested 
widely against various cases of exciton excitations in 
semiconductors,^^ phase diagrams with competing orders 
in transition metal compounds, and iron supercon- 
ductors. ^^^^ A specific combination of the local density 
approximation (LDA) for the first step and DMFT for 
the final step has also been widely applied. They have 
been favorably compared with available experimental re- 
sults. 

However, the accuracy is not clear for more complex 
compounds and for more strongly correlated systems. In- 
deed applications to the organic conductors remain a 
grand challenge. The ab initio low-energy models have 
recently been derived for the a^-ET compounds along the 
line of the above first and second parts of the three-stage 
scheme in the present terminology of MACE.^^ In fact, 
they have shown that the ratio of the typical correlation 
strength (local screened Coulomb interaction), /7, to the 



typical electron transfer t is as large as 10 for the down- 
folded single band models. Furthermore, we have found 
that the ratio of the typical off-site Coulomb interaction 
V to U is as large as 1/3. Elucidating electronic struc- 
tures under such a strong correlation is a challenge at the 
forefront of research to develop efficient and accurate nu- 
merical algorithms. 

The purpose of this study is to make a further step 
to the third procedure following the spirit of MACE; 
we solve the ab initio effective Hamiltonian of real hc- 
ET compounds by using an accurate solver based on a 
recently developed and improved mVMC method with 
many variational parameters'^ to see whether the present 
general framework combined with the mVMC solver of- 
fers an accurate framework for the complex and strongly 
correlated organic conductors. For this purpose, we study 
A>:-NCS, as a typical compound close to the Mott transi- 
tion. 

This compound is in fact barely metallic with an en- 
hancement of the antiferromagnetic correlations revealed 
by the nuclear relaxation rate Ti and located close to 
the antiferromagnetic Mott insulating phase. ^ In fact, 
above 90 K,^^ the resistivity shows insulating-like in- 
crease with decreasing temperatures and has an inflec- 
tion point around T=55 K (ref. 36) coinciding with the 
peak of l/(TiT) with a crossover to a metallic behavior 
below it. It eventually becomes superconducting below 
around 10 K. 

In this study, we assume normal states (not the super- 
conductor) for the candidate of the metal, while leave 
it arbitrary for the insulator to study competitions be- 
tween metals, antiferromagnetic or charge ordered insu- 
lators as well as the Mott insulator without symmetry 
breakings. We examine the phase boundary between the 
metal and the Mott insulator in a parameter space by 
taking the relative electron correlation amplitude as a 
parameter beyond the ab initio value. In practice, we 
draw a phase diagram as a function of a parameter A 
that monitors the relative strength of the effective inter- 
action to the electron transfer by uniformly scaling the 
interaction strength, where the ab initio value is given by 
A = 1. It gives us an idea about the relative location of 
the real material to this phase boundary and hence the 
relevance of the Mott physics. As we will show in §3, off- 
site Coulomb and exchange interactions largely stabilize 
a paramagnetic metal in the region of strong correlation 
{U/t > 6). 

The present work is, to our knowledge, the flrst ab 
initio attempt to estimate the Mott transition and its 
neighboring phases in organic conductors containing a 
large number of atoms beyond 100 with four complex 
ET molecules in a unit cell. Our results indicate that 
A^-NCS is indeed near the Mott transition within the ac- 
curacy of 20%. It also shows that antiferromagnetic in- 
sulating state exists even in the metallic phase near the 
Mott transition, as a metastable excited phase, whose 
energy is typically about 1 meV ~ 10 K higher than 
that of the paramagnetic phase. This is also consistent 
with the above experimental results of the crossover be- 
tween high-temperature insulating and low-temperature 
paramagnetic metallic phases. Although the results show 
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overall agreement with the experiments, a closer look of 
the ground state obtained with the realistic ab initio pa- 
rameter (i.e., A = 1) becomes an antiferromagnetic Mott 
insulator in contrast to the metallic phase in the exper- 
imental indications. We then discuss possible origins of 
this discrepancy. 

This paper is organized as follows: We briefly summa- 
rize in §2 the method of our calculation with ab initio 
downfolding of the low-energy model for /^-NCS and the 
basic framework of the mVMC method we employed. In 
§3 calculated results are presented. The summary and 
discussions are given in §4. 

2. Method 

2.1 Derivation of low- energy effective model 

Here, we describe a derivation of low-energy effective 
models for the present system. The scheme is based on 
first principles calculations and an application of the first 
two stages of the three-stage scheme. The basis of the 
Hamiltonian is the Wannier function associated with an- 
tibonding states of the highest occupied molecular or- 
bitals (HOMOs) of two ET molecules that form a dimer. 
In the present paper, we restrict our consideration within 
the single-band models, where we derive the model con- 
taining only the degrees of freedom for the antibond- 
ing band while the bonding band is traced out in the 
downfolding. Possible dynamical effects of the bonding 
degrees of freedom remain as future issues. The explicit 
form of this Hamiltonian is given in the form of the two- 
dimensional (2D) single-band extended Hubbard model 
as 

a i^j (jp i,j 

'^l^Yl {4a(^]p^iP^j^ + «L4p%P%>t) , (1) 

where aj^ (a^cr) is a creation (annihilation) operator of an 
electron with spin a in the Wannier orbital localized at 
the ith BEDT-TTF dimer. The tij parameters are given 

by 



ti 



il'^Ksl'^i) 



(2) 



with \(f)^) = al\0) and Hks being the Kohn-Sham Hamil- 
tonian representing an effective one-body potential. The 
Vij and Jij parameters are screened Coulomb and ex- 
change integrals in the Wannier-orbital basis, respec- 
tively, expressed as 

Vij = {Mj\W\<j>i<Pj) 



dvdv'4>nr)M^)W{r,r')rAr')4>j{r') (3) 



and 



Jij {(l)i(l)j\W\(l)j(l)i) 



drdrV*(r)(/>,(r)W^(r,rO(/>*(rO</>^(rO (4) 



with W{r^r') being a 2D screened Coulomb interaction 
in the low- frequency limit. Note that the Hamiltonian 



given in eq. (1) can be rewritten as 

cr i/j i 



4e 



(5) 



where the onsite Hubbard parameter U is given by Vu^ 
^io- = ^Icr^io- and Ui = + n^|. The spin operator Si is 
defined as Si = {Sf^S^.S^), where Sf /iS^ = (aj^a^ ± 

al^ai^)/2 and 6'f = (n^^ - nii)/2. 

The derivation of W for the purely 2D system fol- 
lows ref. 37, where a new framework of the constrained 
random-phase approximation (cRPA) was developed for 
the purpose to derive effective interactions of models de- 
fined in lower spatial dimensions. This new scheme is 
suitable for quasi-low-dimensional materials such as the 
present system. The cRPA method is originally formu- 
lated in the RPA framework with the constraint for the 
band degree of freedom to eliminate only the degrees 
of freedom far from the Fermi level in energy. This is 
called the band downfolding. In the proposed scheme of 
the supplementary downfolding,^'' however, the concept 
of the constraint is additionally relaxed to include the 
screening by the polarization in the other layers/chains 
even within the target bands. This is formulated in the 
real space representation and eliminates the degrees of 
freedom away from the target layer/chain, which results 
in the low-dimensional model for the target layer/chain. 
We call it the dimensional downfolding. 

Practically, the band+dimensional downfolding is per- 
formed in two steps: We first perform the band downfold- 
ing to derive the 3D model for small number of bands 
near the Fermi level. This is followed by the dimen- 
sional downfolding in the second step.^^ With this idea, 
we can naturally derive the low-energy model in any di- 
mensions. In the present case, we use it for the derivation 
of a 2D model for a^-SCN. 

2.2 Multi- variable variational Monte Carlo method 

To investigate ground-state properties of the low- 
energy 2D model, we employ a multi- variable variational 
Monte Carlo method (mVMC) combined with quantum- 
number projection and multi- variable optimization.^^' 
The variational wave function |?/;) is defined as 

1^) = ^J^d-h^G>C^=°|0pair), (6) 

where Vq^ ^d-h' Gutzwiller factor, the 

doublon-holon correlation factor, ^^'^^ and the Jastrow 
factor, respectively. These factors are defined as 



Vg = exp 



-g Y rii^riii 



d-h 



exp 



>- m=0^=l,2 



(7) 
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Vj = exp 



'it 



(9) 



where ^[riy ^ind Vij are variational parameters. Here, 
^i{rn) ^ many-body operator which is diagonal in the 
real-space representations. When a doublon (holon) ex- 
ists at the i-th site and m holons (doublons) surround at 
the ^-th nearest neighbor, ^^^^^ gives 1. Otherwise, <ff^^^ 
gives 0. In the present study, we take Vij = vivij) = 
v{—rij)^ where r^j = — Vj is the relative displacement 
between the sites i and j. The spin quantum- number pro- 
jection operator C^^^ restores the SU (2) spin-rotational 
symmetry with the total spin 5* = 0.^^'^^ 

For the one-body part |(/)pair), we employ a generalized 
pairing wave function defined as 



N/2 



Impair) 



|o), 



(10) 



where fij and are variational parameters and the 
total number of the sites, respectively. The number 
of electrons is denoted by N. Throughout this paper, 
we consider the half-filled case TV = TVg. The paring 
wave function given in eq. (10) can flexibly describe 
paramagnetic metals/insulators, antiferromagnetic insu- 
lators/metals, charge-ordered insulators/metals, super- 
conducting phases as well as phases with strong quan- 
tum fluctuations showing developed spin and/or charge 
correlations decaying with any power laws as functions 
of distance. In fact, the mVMC method based on the 
variational function in eq. (6) can describe paramagnetic 
metals with developed spin correlations in hole-doped 
Hubbard models on square lattices. 

In principle, it is better to allow variational wavefunc- 
tions as much as flexible without imposing any constraint 
on fij. However, if we do not impose the constraint, the 
possible variational parameters fij increase in proportion 
to N^, which increases the computation time for the op- 
timization of the variational parameters enormously for 
large system sizes. To reduce the computational cost, we 
restrict fij to have a igx x ^sy sublattice structure as 
fij = fa{j){ri — Tj), where a{j) is a sublattice index 
at site j. In the study of the Mott transition in §3.2.2, 
we take igx = ^sy = 2 (see Fig. 1). This assumption is 
based on the observation of the staggered magnetization 
with ordering vector of (tt, tt) in the frustrated Hubbard 
model at t^t < 0.6,^'^ and the spin-^ quantum Heisen- 
berg antiferromagnets corresponding to U ^ +oo for 
J2/J1 < (0.6)^.^^ These appear to be relevant to the 
present parameter region as we will see in the compari- 
son with the exact diagonalization as well. For the study 
of the charge order in §3.3 (as we describe in detail there), 
we take more general isx = 6 and isy = 2 to allow the 
three-sublattice as well as two-sublattice structures (See 
Fig. 2). 

In the following calculations, we take = L x L sites 
with periodic boundary conditions. All the variational 
parameters are simultaneously optimized by using the 
stochastic reconflguration method. 




Fig. 1. (Color online) Schematic illustration of a 4 x 4 super cell. 
The labels 1-4 denote the sublattice indices for the 2x2 sublat- 
tice structure for fij. The arrows represent the staggered mag- 
netization in the AFI (see Fig. 5). 




Fig. 2. (Color online) (a) Schematic of the three-fold charge- 
ordered (CO) phase. The colored (shaded) and white sites denote 
charge-rich and charge-poor sites, respectively. The labels denote 
the sublattice indices for the 6x2 sublattice structure for fij. 



3. Results 

In this section, we show our computed results including 
ab initio band calculations, derived transfer and interac- 
tion parameters in the low-energy 2D model for /^-NCS, 
and analyses for this model with the mVMC method. 
Our calculation has been performed by using the experi- 
mental lattice structure of /^-NCS taken from the neutron 
diffraction data at 15 K by Schultz et al.^ 

3.1 Low- energy model 

The present ah initio calculations were performed 
with an electronic-structure code based on plane- 
wave basis set, Tokyo Ah initio Program Package.^^ 
Density-functional calculations with the generalized gra- 
dient approximation (GGA) using the Perdew-Burke- 
Ernzerhof parameterization^^ were performed with the 
Troullier-Martins norm-conserving pseudopotentials^^ in 
the Kleinman-Bylander representation.^^ The cutoff en- 
ergies in wavefunctions and charge densities were set to 
36 Ry and 144 Ry, respectively. We employed a 5 x 5 x 5 /c- 
point sampling for the Brillouin-zone integral. The elec- 
tronic structure with the cutoff energy of 36 Ry was com- 
pared with the higher cutoff of 49 Ry. We conflrmed that 
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the 36-Ry band dispersion is almost identical to the 49- 
Ry one. We also confirmed that the wavefunctions of the 
BEDT-TTF molecules in the low-energy, being essential 
in the polarization calculation, are well converged for 
the present cutoff. The construction of the maximally- 
localized Wannier functions follows ref. 51. The polar- 
ization function was expanded in plane waves with an 
energy cutoff of 5 Ry and the total number of bands 
considered in the polarization calculation was set to 750, 
where the numbers of occupied, partially-occupied, and 
unoccupied bands are 232, 2, and 416, respectively. This 
condition corresponds to considering excitations up to 
~21 eV above the Fermi level. The Brillouin-zone in- 
tegral on wavevectors was evaluated by the generalized 
tetrahedron method. The additional terms in the long- 
wavelength polarization function due to nonlocal terms 
in the pseudopotentials were explicitly considered follow- 
ing ref. 53. In the evaluation of the Wannier matrix ele- 
ments, Vij and Jij, the singularity in the Coulomb inter- 
action at the long-wave-length limit was treated in the 
manner described in ref. 53. We confirmed that these 
conditions give well converged results. 

Figure 3(a) shows our calculated GGA band struc- 
ture of A^-NCS. (Blue) dotted lines are a tight-binding 
band with three transfers up to the third nearest neigh- 
bors, listed in Fig. 4(a), derived with matrix elements of 
the Kohn-Sham Hamiltonians in the maximally- localized 
Wannier-orbital basis. The resulting bandwidth of the 
target band for the effective model is 0.56 eV. We note 
that interlayer transfers are considerably small as ~ 0.1 
meV compared to intralayer transfers ~ 65 meV and thus 
the present system is regarded as a typical quasi-2D sys- 
tem. Figures 3(b) and (c) are visuahzation of our calcu- 
lated Wannier function in the side and top views, respec- 
tively. From the figures, we see that the Wannier orbital 
is the anti-bonding state of HOMOs of two ET molecules 
and confined in the layers. 

Figure 3(d) shows the calculated cRPA interactions 
plotted as a function of distance between centers of 
the Wannier orbit als. (Green) filled circles represent the 
conventional cRPA result based on "band downfolding" 
scheme. This cRPA interaction exhibits a power-law de- 
cay (dashed hues) 



at long distances with the dielectric constant e = 5.0 and 
the unit parameter s = 14.40 eV-A. 

On the other hand, (black) triangles in Fig. 3(d) de- 
scribe the result obtained from the present "band + di- 
mensional downfolding" scheme. After this dimensional 
downfolding, the effective interaction becomes qualita- 
tively different from the long-ranged form and is reduced 
to a short-ranged interaction. In fact, it fits well with the 
Yukawa type form (solid lines) as 

er 

where a = 16. 4 A is the interlayer distance between the 
ET layers. This qualitative change into the short-ranged 
interaction comes from the screening by the gapless po- 



larization channel of other metallic layers, which enters 
in the dimensional downfolding process. To explicitly 
distinguish, hereafter, we refer to the former yielding 
eq. (11) as 3D-cRPA and to the latter with eq. (12) as 
2D-cRPA. For comparison, the figure includes the bare 
((blue) squares) and full-RPA ((red) crosses) results as 
well. The dotted line is the bare Coulomb interaction 
decaying as s/r. 

On the basis of this exponential dependence in the 2D- 
cRPA, it is justified to employ the Hubbard- type model 
with only the short-ranged interaction. Practically, Vij is 
considered up to the third nearest neighbors in this pa- 
per, because Vij is negligible beyond this range. We also 
note that the interlayer screening affects even the onsite 
interaction, reducing U=OM eV for 3D-cRPA to [7=0.64 
eV for 2D-cRPA by ~ 25 %. For reference, we note that 
the full-RPA U is 0.19 eV and thus the intralayer screen- 
ing reduces the effective Hubbard U further by 70 %. 

In the calculation of the interlayer screening for the 
dimensional downfolding, we consider stacked supercells, 
each of which contains one target layer/chain. Electrons 
on this target layer/chain are screened by those on other 
layers. The super cell is employed just for the technical 
reason of the calculation and the supercell size should 
be extrapolated to the infinity afterwards. Note that the 
size of the superlattice A^l is nothing but the number of 
the sampling-/c points along the a* axis perpendicular to 
the layer. In the present case, we consider the case where 
the total system consists of the supercell containing up 
to five stacking layers (one target layer and up to four 
screening layers), where we sample five /c-points along the 
interlayer a* axis for the case of the maximum supercell 
size. In Fig. 3(e), we extrapolate Vij to the thermody- 
namic limit of A'l = +CO with an exponential function. 
The obtained Vij in the thermodynamic limit are sum- 
marized in Fig. 4(b). The estimated [//t is as large as 9.6, 
apparently supporting that /^-NCS belongs to a material 
with strongly correlated electrons. 

Screened (direct) exchange interactions Jij are also 
calculated with 2D-cRPA. The results are shown 
in Fig. 4(c). The exchange interactions are hardly 
screened^^'^^ and have no significant system-size depen- 
dence. One might think that Jij are negligibly small com- 
pared to the Coulomb interactions Vij. However, we will 
see that the exchange interactions have a discernible ef- 
fect on the critical interaction ratio of the metal-insulator 
transition as demonstrated in Sec. 3.2. For the nearest 
neighbor pairs in the strong coupling limit, the kinetic ex- 
change is estimated by —2\tij\'^/{U — Vij) as —20.0 meV. 
This value is comparable to the nearest-neighbor direct 
exchange interaction Jij = 6.5 meV. Indeed, Jij and Vij^ 
which are absent in the simple Hubbard model, play a 
substantial role in quantitatively determining electronic 
properties of a>:-NCS as demonstrated in the next section. 

3.2 Ground- state properties of the low- energy model 

We now present the ground-state properties of our de- 
rived 2D low-energy model obtained by using the mVMC 
method after extrapolating finite-size results to the ther- 
modynamic limit. 
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Fig. 3. (Color online) (a) Calculated GGA band structures (red line) of /t-(BEDT-TTF)2Cu(NCS)2. The crystal structure contains 
alternating layers (parallel to the be plane) of BEDT-TTF donor molecules and polymeric Cu(NCS)^ anions. Band dispersions are 
plotted along the high-symmetry points in the be plane, where F = (0, 0, 0), Y = (0, 6*/2, 0), Z = (0, 0, c*/2), and M = (0, 6*/2, c*/2). 
Note that the a axis is interlayer axis. The zero of energy is the Fermi level. The (blue) dotted dispersions are obtained by the three 
transfer parameters listed in Fig. 4(a). The panels (b) and (c) display isosurface contours of maximally localized Wannier functions for 
the target band, where (b) and (c) show the side and top views, respectively. The amplitudes of the contour surface are 0.02 [light grey 
(yellow)] and —0.02 [dark grey (blue)] in the atomic unit, (d) Calculated screened Coulomb interactions of ^c-(BEDT-TTF)2Cu(NCS)2 
as a function of the distance between the centers of maximally localized Wannier orbit als. The (blue) squares, (green) circles, (black) 
triangles, and (red) crosses represent the bare, 3D-cRPA, 2D-cRPA, and full-RPA interactions, respectively. The dotted, dashed, and 
solid curves denote s/r, s/er, and s exp(— r/f7)/er, respectively, where a decay constant e = 5.0 was determined by the fitting to the 
3D-cRPA data. Also, s is a unit parameter of 14.40 eV-Aand a is the characteristic screening length of the 2D-cRPA interaction, 
corresponding to the interlayer distance of 16.4 A. (e) Convergence of interaction parameters Vij up to the third neighbors as functions 
of number of screening layers n contained in a unit cell. Vij are extrapolated to n — oo. Solid curves are the fitted exponential 
functions. Here, V, V' ^ and V" represent the first, second and third neighbor interactions, respectively. 



(a) Transfer integral (b) Screened Coulomb (c) Screened exchange 




Fig. 4. (Color online) Schematic illustration of 2D BEDT-TTF layer, where a BEDT-TTF molecule is described as an ellipsoid and 
its dimer is written as a circle. In the dimer limit, the system forms anisotropic triangular lattice. The derived parameters for the 
2D extended Hubbard model in eq. (1) are shown for transfer integrals (a) and screened Coulomb (b) and screened exchange (c) 
interactions up to the third neighbors. 
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Fig. 5. (Color online) Calculated Xu-^v,J ground-state phase di- 
agram of the model given in eq. (1). PMM, AFI and CO de- 
note a paramagnetic metal, an antiferromagnetic insulator and 
a charge-ordered phase, respectively. The transition points on 
the diagonal line and the horizontal line of Xv,J = are those in 
the thermodynamic limit. The transition points on the horizontal 
lines of Xv,j = 1 and 0.7 are determined by the data for L = 12, 
which are expected to be close to the thermodynamic limit. The 
arrow with dotted line between CO and AFI denotes the transi- 
tion line between uniform and three-fold charge-ordered phases 
analytically estimated in the limit of A[/, Xv,j +oo. 

3.2.1 Ground- state phase diagram 

For a comprehensive understanding of low-energy elec- 
tronic structures of the ab initio model, we work out the 
ground-state phase diagram in the parameter space of 
Xu and Ay^j, where U is scaled from the realistic value 
by the factor A[/, while Vij and Jij are scaled by the 
factor Ay, J. The obtained phase diagram is illustrated 
in Fig. 5. In the weakly-correlated region (namely in the 
region of Af/ <C 1 and Xy^j <C 1), the ground state is a 
paramagnetic metal (PMM). With increasing the inter- 
actions along the diagonal line of A=A[/=Ay,j, the sys- 
tem undergoes, at A:^0.782±0.005, a Mott transition into 
an antiferromagnetic insulator (AFI) with the ordering 
vector at (7r,7r). Thus, the present three-stage approach 
reasonably reproduces the experimental fact that /^-NCS 
is on the verge of the metal-insulator transition. 

The phase diagram tells us that the off-site interactions 
Vij and Jij largely stabilize the paramagnetic metal; for 
the case of Ay^j = 0, the system turns into the Mott 
insulator at Xu — 0.58 ±0.04, which is 20% smaller than 
that for the diagonal line. On the other hand, a charge- 
ordered (CO) phase emerges for Xu ^ Ay^j. In the CO 
phase, the system exhibits a three- fold {rich-poor-poor) 
charge order with ordering vector of (27r/3, 27r/3), which 
reduces the nearest-neighbor Coulomb repulsion energy 
(see an inset of Fig. 5). The nature of the CO phase will 
be discussed in §3.3. 

3.2.2 Mott transition 

In this subsection, we discuss the Mott transition be- 
tween the paramagnetic metal and the antiferromagnetic 



insulator. Let us start with analyses on the diagonal line 
of Xu = Ay, J. To identify quantum phase transitions, we 
calculate the doublon density 

^ i 

the momentum distribution n{k) 

"(fc) ^ ^ E<4c,<.)e-"^-(--'-^), (14) 
and the spin structure factor 

^(9) ^ ^J2(^^■S,)e-"'<'■'-^^\ (15) 

In Fig. 6, we compare calculated results obtained by 
the unrestricted Hartree-Fock (UHF), and mVMC cal- 
culations with those by exact-diagonalization methods 
for L = 4, to roughly get insight into accuracies of the 
mVMC and UHF methods. Hereafter we measure sys- 
tem sizes in units of the unit cell consisting of a square 
constructed from the nearest-neighbor bonds, in which 
the linear dimension is denoted by L. In all of the three 
methods, the results show qualitatively similar behav- 
iors, where by turning on A, I) decreases continuously 
from 0.25 and suddenly exhibits a discontinuous decrease 
at Ac [panel (a)]. At the same time, S{Q)/Ns shows 
a jump at A = Ac [(b)]. For A > Ac, a clear single 
peak emerges at Q = (7r,7r) in S{q)/Ns as shown in 
Fig. 6(c). We found no other anomaly in S{q) and D. 
These numerical results suggest that the system is under 
the influence of the electron correlations and undergoes 
a first-order transition from the paramagnetic metal to 
the antiferromagnetic insulator at Ac- For the exact di- 
agonalization, we obtain Ac ~ 1.49 for L = 4. On the 
other hand, the mVMC result indicates Ac = 1.15-1.2, 
which is about 20% smaller than the exact value. De- 
spite the underestimate of the critical interaction ratio, 
the mVMC well reproduces S{Q)^ and n{k) in the 
AFI and PMM phases [see also Fig.6(d)]. It should be 
mentioned that a UHF calculation gives Ac = 0.6-0.65, 
which is over 50% smaller than the exact value. To sum- 
marize, the present small cluster result implies that the 
mVMC solver tends to underestimate Ac by as much as 
20% compared to the exact-diagonalization solver. The 
comparison between the mVMC and PIRG results for 
the Hubbard model on the square lattice with the next- 
nearest-neighbor transfers has shown the underestimate 
of Ac for the metal-insulator transition by the mVMC as 
much as 10 % in the thermodynamic limit. We further 
discuss the implication of this discrepancy in §4. 

Figure 7 shows system-size dependences of D [panel 
(a)] and S{Q)/Ns [(b)] up to L=12, obtained by the 
mVMC calculation. We found that these physical quanti- 
ties are well converged for L > 6. For all the system sizes, 
the Mott transition is characterized by a level crossing 
between the AFI and PMM states, and a finite jump in 
the properties, indicating their first-order nature of the 
transitions. 

We next consider the ground-state properties in the 
thermodynamic (bulk) limit. To this end we extrapo- 
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Fig. 6. (Color online) Calculated doublon density D and spin structure factor S(Q) for L = 4. Here Q denotes (7r,7r). For details, see 
text, (c) Spin structure factor S(q) and (d) momentum distribution n(k) at A = 1.0 and 1.5. 
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Fig. 7. (Color online) System-size dependences of (a) doublon 
density D, and (b) spin structure factor S(Q). Here Q denotes 
(7r,7r). (c) D, mg, and mg extrapolated to the thermodynamic 
limit. All the data are for Xjj = Xv,j which corresponds to the 
diagonal line in Fig. 5. 



late the energies per site E/N^ of the AFI and PMM 
states around Ac- For the AFI state, we employ the scal- 
ing form AE/Nq (x given by the spin- wave theory 
for the two-dimensional quantum Heisenberg antiferro- 
magnets^^ with AE being the finite-size correction of 
the total energy. For the PMM state, we extrapolate the 
energy by following AE/N^ oc The top and mid- 

dle panels in Fig. 8(a) illustrate the procedure of the size 
extrapolation for PMM and AFI, respectively. The criti- 
cal interaction ratio was estimated as Ac = 0.782 ± 0.005 
[the bottom in Fig. 8(a)] as the crossing point of the two 
energies extrapolated to the thrmodynamic limit. 

After determining Ac, we calculated the ground-state 
physical quantities in the thermodynamic limit. The dou- 
blon density D in the thermodynamic limit for each lo- 
cally stable state was estimated with the scaling form 



D{L = oo) — D{L) (X L~^. We also calculated the stag- 
gered magnetization mg by extrapolating S{Q)/Ns to 
the thermodynamic limit as mg — S{Q^L)/Ns oc L~^. 
This scaling form for mg is suggested by the spin-wave 
theory for the two-dimensional quantum Heisenberg an- 
tiferromagnets.^^ The extrapolated I), mg, and mg are 
shown in Fig. 7(c). Finite jumps in D and mg in the 
thermodynamic limit indicates the first-order nature of 
the Mott transition, which is consistent with the experi- 
mental behavior observed in the pressure-controlled Mott 
transition of tv-Cl.^^ The obtained ordered moment is 
rus ~ 0.22 near the metal-insulator transition. This value 
is much smaller than that for the spin- 1/2 Heisenberg 
model on the square lattice (- 0.3),^^ while it IS com- 
parable to that in the antiferromagnetic phase near the 
Mott transition in the geometrically frustrated lattice. 

In Fig. 9(a), we plot the momentum-resolved spin 
structure factor S{q) calculated for L = 12 near the 
metal-insulator transition. There is no essential differ- 
ence between the result for A = and that for A = 0.75, 
indicating the absence of the antiferromagnetic spin fluc- 
tuations near the metal-insulator transition. This is also 
seen more directly in the absence of system-size depen- 
dence for S{Q) [Fig. 9(b)]. The antiferromagnetic cor- 
relation length ^AF was estimated from a fitting of the 
Ornstein-Zernike form 

S{Q) 



S{q) 



(16) 



i + aF(9-Q)' 

to have ^af/<^ = 0.42 ±0.02 with a being the lattice spac- 
ing in terms of the square lattice [see Fig. 9(c)]. This is a 
strong indication that the antiferromagnetic correlation 
does not develop even near the metal-insulator transi- 
tion, which is consistent with the experimental observa- 
tion of K-NCS at low T; i.e., T < 55 K.^^ We note that 
the enhancement of the antiferromagnetic spin correla- 
tions above 55 K in /^-NCS may be due to proximity 
effects of the Mott transition. For instance, at A = 0.75 
(about 4% below Ac, which may be a typical plausible 
value for the real compound), the energy difference per 
site between AFI and PMM is on the order of 1 meV 
10 K, and thus it is likely that the antiferromag- 
netic metastable state partially contributes to the finite- 
T manifold and enhances the antiferromagnetic correla- 
tions. Further studies for finite-T properties are desirable 
for this topic. 

A recent mVMC study indicates that the charge exci- 
tation gap is partially formed in the metallic phase as a 
precursor of the Mott gap for the square-lattice Hubbard 
model with next-nearest hopping, where an "arc-like" 
Fermi surface is observed near the Mott transition, i.e., 
Uc/t = 3.3 at t^t = —0.3. To make a comparison with 
this result, we plot in Fig. 10 the momentum distribu- 
tion n{k) and its gradient |Vn(/c)| for A=0.75 [panel (a)] 
and 0.8 [(b)]. The Fermi surface of A = is denoted by 
dotted lines in the contour plot. The electron correlations 
smear the jumps in n{k) and this effect is more significant 
around (tt, 0) than around (7r/2,7r/2) [Fig. 10(c)]. We re- 
mark that the present result has no "arc-like" structure 
in the \\/n{k)\ plot near the Mott transition in contrast 
to the previous study. This may be due to the large crit- 



10 J. Phys. Soc. Jpn. 



Full Paper 



H. Shinaoka et al. 



ical interaction strength, i.e, Uc — 7.5, and the resultant 
strong first-order nature. ^'^ 

To discuss effects of the off-site interaction Vij and 
Jij , we compare the critical interaction ratio Ac obtained 
with/ without these parameters. By switching off Vij and 
Jij (i.e., the line along Xy j=0 in Fig. 5), we obtained 
Ac=0.58±0.04 [Fig. 8(b)]. By switching off only Jij (i.e., 
the Xu=^v,j line but Jij=0), we obtained Ac=0.74±0.01 
[Fig. 8(c)]. Since Ac is 0.782±0.005 for the model having 
Vij and Jij, we have an appreciable enhancement of Ac 
with Vij and a minor modification by Jij. The mech- 
anism of the increase in Ac by Vij is as follows: With 
introducing Vij , the creation energy for a doublon-holon 
pair in the Mott insulator is reduced from /7=0.64 eV 
to U — =0.45-0.48 eV. Furthermore, Vij also reduces 
the local moment miocai = \/l — 2D/2 by increasing the 
doublon density. The introduction of J^j, on the other 
hand, favors a ferromagnetic correlation on the nearest- 
neighboring bonds, thus destabilizes the antiferromag- 
netic solution. 

3.3 Charge- ordered phase 

In general, off-site Coulomb interactions tend to sta- 
bilize charge-ordered states. This effect is not fully con- 
sidered in the previous section, in which the analyses are 
limited to 2 x 2 orders; e.g., at half filling, the nearest- 
neighbor repulsion stabilizes three- fold charge orders. In 
this section, we examine the possibility of a charge or- 
dering induced by Vij in more extended parameter space 
with a more generalized variational wave function. 

For calculations along the lines of Xy^j = 1 and 
^v,j = 0-7 in the phase diagram of Fig. 5, we extend 
the pairing wave function Impair) by allowing fij to have 
a 6 X 2 sublattice structure (see Fig. 2). This allows us to 
search magnetic and charge orders with ordering vectors 
of (tt, tt) and (27r/3, 27r/3) on equal footing. To identify 
charge orders, we calculate the charge structure factor 

^(9) = JTs ^^""^ ~ " l)^""'-^'''"'''^- (17) 

The L=6 and L=12 results calculated for Xy^j = 1.0 
and 0.7 are shown in Fig. 11. By decreasing Xu with 
Xy^j fixed, the system turns into the CO from the PMM. 
At Ay, J = 1.0 and 0.7, this transition is found to be of 
the first order as seen in the panel (a). The CO phase 
is characterized by a three- fold {rich-poor-poor) order as 
seen in the {rii) profile [the bottom of the panel (a)] cal- 
culated at Xu = 0.0 and Xy^j = 1.0. The left side of 
the panel (b) displays our calculated n(/c), |Vn(/c)|, and 
N{q) at Xu=0.0 and Ay,j=1.0 for L = 12. We found 
sharp peaks in N{q) at (27r/3, 27r/3) and (47r/3, 47r/3), 
due to the three- fold CO order. A UHF calculation indi- 
cates that the CO state has a small charge gap of about 
0.2 eV associated with a two- fold bond order at A[/=0.0 
and Ay,j=1.0 for I/=24 (not shown). However, we could 
not conclude whether or not the CO phase is insulating 
in the mVMC calculations [see Fig. 11(b)], because the 
system sizes in the present study are not sufficiently large 
for the size extrapolation of the charge gap. 

A remarkable observation here is that the peaks in 



N{q) exist even in the PMM phase as seen in the right 
side of Fig. 11(b). We estimate a correlation length of the 
charge fluctuations ^c in terms of the Ornstein-Zernike 
function 

with gmax = (2/37r, 2/37r). The result is shown in 
Fig. 11(c): We obtained ^c/a=1.8±0.6 and 1.5±0.3 for 
Af/=0.65 and Xu=0.75 (PMM), respectively, at Ay,j = 
1.0. These correlation lengths are appreciably longer than 
those of the antiferromagnetic correlations; we obtained 
aF/a=0.33±0.01 and 0.34±0.01 for At/=0.65 and 0.75, 
respectively. These indicate that the charge fluctuations 
are more dominant in the PMM phase than the antifer- 
romagnetic correlations. More interestingly, in the PMM 
phase at Xy^j = 1.0 and Af/ = 0.75, the Fermi surface be- 
comes smeared around (±7r, 0) and (0, ±7r), and as a con- 
sequence, "Fermi arcs" appear around (±7r/2, ±7r/2) [see 
the right panels of Fig. 11(b)]. This is in sharp contrast 
to Figs. 10(a), where we did not see the "arc-like" struc- 
ture in \\/n{k)\. This indicates that the emergence of 
the arc-like structure is ascribed to the enhanced charge 
fluctuations. Recently, Fermi- arc like behavior has been 
observed in a metallic phase adjacent to a CO phase for a 
metallic layered nickelate Eu2-ajSra;Ni04.^^ The forma- 
tion of a Fermi arc might be characteristic to metallic 
states with strong charge fluctuations. 

The results are summarized in Fig. 5 as the Xu-Xv,j 
phase diagram. The PMM is sandwiched by CO and AFI 
and 40% reduction of U in PMM realized in Ay,j=l 
causes the transition to CO. The CO transition line and 
AF transition line merge in the strong correlation limit, 
i.e., Ac/, Ay, J +00. The critical interaction ratio in this 
limit is given by Xu/Xv,j = y^y^/2-y^^ — ^•'^^ based on 
a simple estimate of the ground-state energy (see Ap- 
pendix). 

4. Summary and discussion 

We have performed mVMC calculations for an ab ini- 
tio low-energy effective model of /^-NCS. The ground 
state of this compound has turned out to be close to 
the Mott transition. Within 20% change in the ab ini- 
tio interaction parameters, it undergoes a transition be- 
tween an antiferromagnetic insulator and a metal. The 
real compound is known to be scarcely metallic at low 
temperatures while semiconducting above 90K with en- 
hanced antiferromagnetic correlations. The calculated re- 
sult reproduces these basic experimental results. In this 
first challenge of the ab initio calculation of a>:-NCS, the 
present result has proven a good accuracy of the three- 
stage ab initio scheme even for strongly correlated and 
complex organic compounds. 

However, strictly speaking, the real compound of tz- 
NCS becomes metallic at low temperatures implying that 
the present result shows approximately 20% underesti- 
mate of the transition parameter in terms of the interac- 
tion strength. A possible origin of this discrepancy may 
be the overestimate of the antiferromagnetic region by 
the mVMC method in the strong correlation regime. Ac- 
tually, the benchmark of the metal-insulator transition 
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Fig. 8. (Color online) Extrapolation of energy E/Ns to the thermodynamic limit around the Mott transition along (a) diagonal line of 
Xjj = Xv,J in the phase diagram (Fig. 5), (b) horizontal line of Xv,j = 0, and (c) diagonal line of Xjj = Xv,j with Jij = 0. The energy 
E is given in the unit of meV. The top and middle panels illustrate the size extrapolations of the PMM and AFI, respectively. The 
bottom panels show energy comparisons in the thermodynamic limit for the PMM and AFI phases. 
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20.0 



Fig. 9. (Color online) (a) Momentum-resolved S(q) at A 
0,0.75, and 0.8 for L = 12. (b) A dependence of S{Q). 
Ornstein-Zernike type fit of S{q) at A = 0.75 (see text). 



(c) 



for the Hubbard model on the square lattice with the 
next-nearest-neighbor transfer t' = 0.2 ~ 0.3 has shown 
nearly 10% overestimate of the antiferromagnetic insulat- 
ing phase as compared to a more precise PIRG result (the 
metal-insulator boundary is around U/t = 3.6 by the 
PIRG estimate against the mVMC result U/t - 3.3). '^'^^ 
Therefore, in more strongly correlated region as in this 
case in comparison to the benchmark, it is likely to over- 
estimate the stability of the antiferromagnetic insulator 
as well in a similar rate; in fact, 20% overestimate has 
been confirmed in the present comparison with the ex- 
act diagonalization for a small cluster (L=4). The de- 
velopment of a more accurate low-energy solver is a fu- 
ture important issue. In the PIRG, it is possible to esti- 
mate the boundary without an explicit bias beyond the 
present mVMC and it is worthwhile to do so, although 
the present parameter with a large frustration and inter- 
action may demand a heavy calculation. 

Another possible origin of the overestimate of the sta- 
bility of the antiferromagnetic insulator is dynamical 
effects of the HOMO bonding band eliminated in the 
present downfolding procedure. If the dynamical polar- 
ization arising from the HOMO bonding band becomes 
important, one has to consider the two-band model as the 
low-energy effective model. Charge fluctuations within 
the dimer of the BEDT-TTF molecule may melt the anti- 
ferromagnetic order and shift the Mott transition bound- 
ary to a higher U value. This may be related to the di- 
electric anomaly recently observed in the experiments 
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Fig. 11. (Color online) (a) Calculated A^((7max), S(Q), and D at Ay^j = 1.0 and 0.7 as functions of Xu- A^(gmax) is the peak value 
of the charge structure factor N(q). The open and filled circles denote data for L = 6 and 12, respectively. For all Xu and Xv,j, we 
found that gmax = (27r/3, 27r/3), and equivalently at (47r/3, 47r/3). The bottom distorted square lattice denotes the charge density (rii) 
calculated at Xu = 0.0 and Xy^j = 1.0 for system size L = 12. (b) Momentum-resolved n(k), \Vn(k)\ and N(q) at Xv,j = 1-0 and 
L = 12. The broken (yellow) lines denote the Fermi surface at Xu = Xv,j = and L = oo. (c) Ornstein-Zernike type fit for N(q) along 
the (0, 0)-(27r/3, 27r/3) line. The system size is L = 12. 
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quantum spin liquid for k-GN as well as the unconven- 
tional Mott transition for ti-Cl by using the present ac- 
curate ah initio method. 

Before closing this paper, we make a brief discussion on 
a recent experimental observation of a relaxer-like dielec- 
tric anomaly for dimer-Mott insulator /^-CN,^^'^^ which 
is known as a spin-liquid material. As its possible origin, 
effects of charge fluctuations within the ET dimers have 
been extensively investigated theoretically.^^' ^^'^^ These 
effects are not taken into account in the present study 
in which the ah initio model is constructed based on the 
dimer basis. However, the present study has revealed that 
charge fluctuations are also enhanced by the inter-dimer 
Coulomb interactions in a wide range of the parameter 
space, especially near the charge-ordered phase, which is 
not very far from the realistic choice of the parameters 
Af/ = Ay = 1. This indicates that the off-site (inter- 
dimer) Coulomb interactions may also play a certain 
role in the dielectric anomaly observed in /^-CN. First- 
principles studies on a^-CN are left for future study. In 
particular, it is of great interest to understand effects 
of inter/intra-dimer charge fluctuations based on an ah 
initio two-band model of /^-CN. 



Fig. 10. (Color online) Momentum distribution n{k) and contour 
plots of |Vn(fc)| at A=0.75 (a) and 0.8 (b), respectively (L = 12). 
Here, n{k) is calculated by using the value interpolated by the 
bi-cubic interpolation of n(k). The broken lines denote the Fermi 
surface at A = and L = oo. (c) A dependence of n(k) along 
symmetry lines. 



around 6K of a^-CN^^'^^ as discussed below. 

The insulating side of this compound is shown to have 
a remarkably suppressed ordered moment (~ 0.22) of the 
antiferromagnetic order. This is clearly ascribed to the 
geometrical frustration effect as well as the proximity of 
metals. A quantum spin liquid is stabilized for more frus- 
trated case as /^-CN, and this compound is also barely 
antiferromagnetic, suggesting that the ordered moment 
for X = Cu[N(CN)2]Cl (a very similar compound but 
just in the insulating side) is likely to have a small or- 
dered moment similar to the present estimate. 

By applying ah initio downfolding scheme to /^-NCS, 
we have shown that the nearest-neighbor and next- 
nearest-neighbor off-site Coulomb interactions V of the 
2D low-energy effective model are unexpectedly large, 
i.e., V/U ~ 1/4. These large off-site Coulomb interac- 
tions stabilize the paramagnetic metal through the exci- 
tonic effect. 

Experimentally, /t:-NCS is known to be superconduct- 
ing at low temperatures. In this paper, we have assumed 
only the normal state in the metallic phase and have not 
examined the possibility of superconductivity in detail to 
focus on the metal-insulator boundary. The phase com- 
petition involving the superconducting state is left for 
future study, because it requires very subtle comparison 
of the stability of the superconducting state. 

It has been proven that the three-stage scheme is pow- 
erful for the ah initio calculation of the organic conduc- 
tors. It is an intriguing issue to study the nature of the 



Acknowledgments 

The authors thank Daisuke Tahara for the use of his 
mVMC code and useful comments. Numerical calcula- 
tion was partly carried out at the Supercomputer Cen- 
ter, Institute for Solid State Physics, Univ. of Tokyo. 
This work was supported by Grant-in-Aid for Scien- 
tific Research (No. 22740215, 22104010, 22340090, and 
23110708), from MEXT, Japan. A part of this research 
has been funded by MEXT HPCI Strategic Programs for 
Innovative Research (SPIRE) and Computational Mate- 
rials Science Initiative (CMSI). 

Appendix: Critical interaction ratio of the 
three-fold charge-order transition in 
the classical limit 

We assume that the ground state is non-magnetic and 
has a three- fold charge order. In the limit of A[/, Xv,j 
+00 (classical limit), the energy per site E/N^ is given 

by 

E/N, (A-1) 
= XuEu + \v,jEv + \v,jEv + \v,jEv", 



(n\ -\- nl -\- nl) + const. 



--)} 



X (A-2) 



where 



U 



(A-3) 



Ev = ^{9n^- {nl + nl+nl)}, (A-4) 
Ev = y {9n2 - {nl + + 4)} , (A-5) 



V" 

Ev" = + 



(A-6) 
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Here, ni, n2 and 77,3 are the site occupancies for the three 
sublattices, and the mean site occupancy is given by n = 
(ni + n2 + ^3)73 = 1. The Cauchy-Schwarz inequahty 
tehs us that ni +77.2 + 77,3 has it minimum 1/3 when ni = 
n2 = ns = 1/3. This indicates that a three- fold charged- 
ordered state becomes more stable than the uniform state 



when Xu 



A 



V 



v 



V" < 0. Therefore, the 



critical interaction ratio of the three-fold charge-order 
transition is given by 



\u/\ 



V,J 



u 

4 



V 



2 



(A-7) 



in the classical limit. 
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